Characterizing the network topology of the energy landscapes of atomic clusters 
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By dividing potential energy landscapes into basins of attractions surrounding minima and link- 
ing those basins that are connected by transition state valleys, a network description of energy 
landscapes naturally arises. These networks are characterized in detail for a series of small Lennard- 
Jones clusters and show behaviour characteristic of small-world and scale-free networks. However, 
unlike many such networks, this topology cannot reflect the rules governing the dynamics of network 
growth, because they are static spatial networks. Instead, the heterogeneity in the networks stems 
from difl'erences in the potential energy of the minima, and hence the hyperareas of their associated 
basins of attraction. The low-energy minima with large basins of attraction act as hubs in the net- 
work. Comparisons to randomized networks with the same degree distribution reveals structuring 
in the networks that reflects their spatial embedding. 
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I. INTRODUCTION 

In recent years, characterizing the energy landscape of 
a system in order to gain a better understanding of its 
behaviour has become an increasingly popular research 
approach,^ with many applications in the fields of protein 
folding,^ clusters and supercooled liquids. '^"^ For exam- 
ple, a common explanation for how proteins are able to 
overcome the Levinthal paradox^ and fold to their native 
state on biologically reasonable time scales is in terms of 
a funnel-like energy landscape. 

Here, we want to take a somewhat different approach 
and focus not on the relationship between the landscape 
and a system's behaviour, but on some of the fundamen- 
tal properties of potential energy landscapes. For exam- 
ple, the number of minima increases exponentially with 
system size^"^ and the size scaling for higher-order saddle 
points has also been obtained. ^'^"^^ It is also known that 
the distribution of minima should be a Gaussian function 
of the potential energy.^'^'^'' 

In this paper we seek to provide fundamental new in- 
sights into the structural organization, and particularly 
the connectivity, of an energy landscape. To achieve this 
we first map the landscape onto a network, and then 
analyse the topology of this network for a series of small 
clusters for which complete networks can be obtained. A 
brief report of some of this work has already appeared. 

Our analysis of these networks is heavily influenced by 
the literature on complex networks. ^^"^^ The recent in- 
terest in this area was in part sparked by the seminal 
paper of Watts and Strogatz,^^ who showed that many 
real-world networks have features that are typical both 
of a lattice (the presence of local order) and of a ran- 
dom graph^^^^^ (short average separations between the 
nodes), and introduced a 'small- world' network model 
that could represent both these features. Subsequently, it 
has been shown that many networks also have a probabil- 
ity distribution for the number of connections to a node 
(in network parlance, the degree distribution) that has 
a power-law tail. Such 'scale-free' networks^^ have been 



found in an impressively diverse range of fields, including 
astrophysics,^^ geophysics,^'' information technology,^^ 
biochemistry,^^'^° ecology"^^ and sociology. "^^ 

The paper is organized as follows. In Section II we ex- 
plain how energy landscapes can be described in terms of 
networks. In Section III we present a detailed character- 
ization of these networks. This analysis initially focusses 
on whether the networks show small-world and scale-free 
behaviour, before going on to look at more subtle features 
of these networks. Then, in Section IV we discuss the im- 
plications of our results for understanding the dynamics 
on complex energy landscapes. 



II. ENERGY LANDSCAPES AS NETWORKS 

To model an energy landscape as a network one must 
first decide on a definition both of a state of the sys- 
tem and how two states are connected. The states and 
their connections will then provide the nodes and edges 
of the network. For systems with continuous degrees of 
freedom, perhaps the most natural way to achieve this 
is through the 'inherent structure' mapping of Stillingcr 
and Weber. ■^'^ In this mapping each point in configura- 
tion space is associated with the minimum (or 'inher- 
ent structure') reached by following a steepest-descent 
path from that point. This mapping divides configu- 
ration space into basins of attraction surrounding each 
minimum on the energy landscape, and is illustrated for 
a model two-dimensional energy landscape in Fig. 1. 

For systems with large numbers of degrees of freedom, 
the energy landscape is an extremely complicated multi- 
dimensional function. The inherent structure approach 
provides a way of dealing with this complexity by break- 
ing the landscape up into more manageable pieces. For 
example, good approximations to the partition function 
of a basin of attraction and to the rates of transitions 
between two basins can be obtained within the harmonic 
approximation (and if necessary with additional anhar- 
monic corrections), thus allowing thermodynamic and 
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FIG. 1: (Colour online), (a) A model two-dimensional energy 
surface, (b) A contour plot of this surface illustrating the 
inherent structure partition of the configuration space into 
basins of attraction surrounding minima. The basin bound- 
aries are represented by thick lines, and the minima and tran- 
sition states by dots, (c) The resulting representation of the 
landscape as a network. 



dynamic properties of a system to be calculated.^ 

The inherent structure approach also provides a nat- 
ural partition for the dynamics of a system. At suffi- 
ciently low temperature, the system will spend most of its 
time oscillating within a basin of attraction surrounding 
a minimum with occasional hops between basins along 
a transition state valley. '^^ This intcrbasin dynamics can 
then be represented as a walk on a network, whose nodes 
correspond to the minima and where edges link those 
minima that are directly connected by a transition state 
(a first-order saddle point). Such a network is illustrated 
for the model surface in Fig. 1, and we will term them 
inherent structure networks. 

In chemistry such networks are often called reaction 
graphs, and have been used to analyse isomerization in 
small molecules and clusters. In these applications, each 
permutational isomer of a structure is usually consid- 
ered separately. For example, in Ref. 34 the reaction 
graph connecting the eight permutational isomers of wa- 
ter dimer is depicted. Such knowledge is for example im- 
portant for understanding the pattern of splittings that 
results from tunnelling between these isomers. 

For an atomic cluster with a single atom type the num- 
ber of permutational isomers of a minimum is 2N\/h 
where h is the order of the point group of that minimum. 
This factor makes consideration of the network with a 
node for each permutational isomer unfeasible for all but 
the smallest sizes. For example, for a 7-atom Lennard- 
Jones (LJ) cluster the four geometrically distinct minima 
give rise to 8904 permutational isomers, but for LJ14 the 
4196 geometrically distinct minima give rise to 6.68 x 10^^ 
permutational isomers. Hence, here we consider the net- 



work where each geometrically-distinct minimum is rep- 
resented by a single node. Indeed, this approach is per- 
fectly reasonable, because one is usually only interested 
in the structure of the system, and not in the permuta- 
tional order. The network we consider here will also be 
unweighted, i.e. we are only interested in whether two 
minima are connected, and not in the weight, however 
defined, of this connection. 

In the first instance, self-connections and multiple 
edges are also excluded, even though there are transi- 
tion states that mediate degenerate rearrangements be- 
tween permutational isomers of the same minimum, and 
there can be more than one geometrically distinct tran- 
sition state connecting two particular minima. As a con- 
sequence of these exclusions the number of edges in our 
networks, so defined, is roughly 30% less than the total 
number of transition states for the larger clusters (See 
Table I). We exclude self-connections, both because we 
are not interested in permutational order, and because 
degenerate rearrangements make no contributions to the 
structural dynamics, and multiple edges because we are 
not so much interested in the number of connections be- 
tween two minima, but just whether they are connected. 
Besides, if we were concerned with the number of connec- 
tions it would not be as simple as counting the number of 
geometrically distinct transition states that connect the 
two minima because of symmetry factors. The number 
of versions of a transition state that connects to a mini- 
mum for a non-degenerate rearrangement is given by the 
factor ft-min/^ts, 1-6. the ratio of the orders of the point 
groups of the minimum and the transition state. ^ 

The properties of these types of networks have ef- 
fectively been studied before, albeit in a directed form 
where the links have been weighted by the rate of tran- 
sition from one minimum to another at a particular 
temperature. ^^"""^ This approach leads to a rate matrix 
from which an exact solution of the interminimum dy- 
namics of the system can be obtained in terms of the 
eigenvalues and eigenvectors of this matrix using a mas- 
ter equation method. However, in all these applications 
the emphasis was on the resulting dynamics and not on 
the topological structure of the networks. Moreover, be- 
cause the focus was the dynamics of a particular transi- 
tion and not the global dynamics, most of the networks 
studied were incomplete. 

We should mention three other studies that have 
sought to characterize the connectivity of configuration 
space in terms of complex networks. ^^""'^ The first was of 
a two-dimensional lattice polymer, where each discrete 
configuration was a node in the network and configura- 
tions were linked to those that could be obtained by the 
application of a single elementary move using a move set 
of corner flips and "crankshaft" rotations. ^^■''^ The sec- 
ond was of RNA considered at the secondary structure 
level, where a node corresponded to those configurations 
that had the same map of contacts between nucleotides 
and nodes were linked to those that could be obtained 
by an elementary change in the contact map.^^ The main 
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TABLE I: Properties of the energy landscapes and inherent structure networks for LJat for iV = 7 — 14. nmin is the number 
of minima; rits is the number of transition states; riodgcs is number of edges in the network; p is the average probability of an 
edge between any two nodes; (k) is the average degree; fcmax is the maximum degree; /avc is the average separation between 
nodes; r is the radius of the network; d is the diameter of the network; C'l and C2 arc different definitions of the clustering 
coefficient; and ru and rs are the assortativity coefficients^^ with respect to the degree and the potential energy of the 
minimum, respectively. The uncertainties in the last two quantities have been calculated using the jack-knife method. 
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diffc;rc;nccs between these networks and our own is in the 
discrete, rather than continuous, nature of the configura- 
tion space and the definition of a state — there is no equiv- 
alent of the inherent structure mapping that groups con- 
figurations into a single state (or basin). We will discuss 
the consequences of these differences later in the paper, 
when we report the properties of the inherent structure 
networks. The third study, which examined the connec- 
tivity of the conformation space of polypeptides,^* is in 
a much more similar spirit to the current work, because 
the space is continuous and connections have a dynami- 
cal basis. Given the niuc;li larger size of the configuration 
space than in the current study, a more coarse-grained 
picture of the landscape is required. Rao and Caflisch 
probe the connections between free energy minima de- 
fined through structural order parameters that assign a 
secondary structure character to each amino acid. The 
sampling of the landscape is not complete, but instead 
based on long molecular dynamics simulations at a tem- 
perature near to the folding transition. The properties 
of these networks are much more similar to the present 
networks. 



III. NETWORK PROPERTIES 

We chose to characterize the inherent structure net- 
works of a series of Lennard-Jones clusters, whose poten- 
tial energy is given by 
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where e is the pair well depth and 2^^^cr is the equi- 
librium pair separation. We concentrate on small sys- 
tems for which we are able to characterize the energy 
landscape completely, because we are interested in the 
global topology of the complete inherent structure net- 
work. We choose to study clusters, rather than a bulk 



system to avoid complications associated with boimdary 
conditions. For clusters we do not have to worry about 
the density dependence of the networks. Moreover, for 
systems that arc small enough for the complete inher- 
ent structure network to be sampled, periodic boundary 
conditions represent a significant perturbing constraint 
on the system, and so do not provide an accurate repre- 
sentation of a bulk system. By contrast, small LJ clusters 
are interesting systems in their own right, whose struc- 
tural, thermodynamic and dynamic properties have been 
well characterized. 

A (near-)complete sampling of the relevant stationary 
points of LJjv up to 7V=13 had previously been obtained 
in a work looking at the scaling behaviour of the numbers 
of stationary points. In addition we generated samples 
of minima and transition states for LJ14. 

Here, we quickly summarize the approach used. It is 
an iterative procedure whereby starting from each min- 
ima in our sample we perform a series of transition state 
searches. We then step off any new transition states along 
the direction of negative curvature to find the connect- 
ing minima. More transition state searches are then per- 
formed from any new minima. This procedure is repeated 
until the number of transition states appears to have con- 
verged with respect to increasing numbers of searches. 
However, to find all the transition states we also have 
to perform searches for second-order saddle points start- 
ing from the transition states and minima. Prom these 
second-order saddle points we then repeatedly perform 
transition state searches. With this additional tool the 
number of transition states increases significantly and 
does appear to converge to the true total, although of 
course there is no proof of this. 

The number of minima and transition states found by 
this procedure for each cluster size are given in Table 
I. The factor that limits the size for which such com- 
plete stationary point samples can be generated is not 
so much the number of stationary points, but the huge 
number of searches that need to be performed to achieve 
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convergence. We were able to map out the complete in- 
herent structure networks for all clusters up to A^=14. 
The network data is available online.^'"' 

There is a theoretical expectation that the number of 
minima and transition states scale with size as^'^'^*^ 



and 



nts 



aNe 



aN 



(2) 



where a and a are constants. These forms are based on 
the assumption of a bulk system (i.e. large A'' and no 
surface) that can be divided up into independent equiv- 
alent subsystems. Hence, one should not expect these 
forms to hold too precisely for small clusters. Expressions 
for the number of higher-order saddle points can also be 
derived. ^^'^^ Based on the Hessian index for which the 
number of stationary points is a maximum, a value of 
a=2 appears appropriate for clusters. 

In network parlance, the number of connections to a 
node is called the degree and denoted by k. Taking the 
above forms for nmin and nts as appropriate and ignoring 
for the moment that not all transition states give rise 
to edges in the networks (because of our exclusion of 
multiple edges and self-connections), one expects that 
the average degree for our networks should follow 



(fc)«^=2a7V=— logn„ 

TT-min OC 



(3) 



i.e. {k) should scale linearly with the size of the clusters, 
or logarithmically with the size of the network. It follows 
that p, the fraction of all possible edges in the network 
that are actually connected is given by 
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Hence the networks should become more sparse as they 
increase in size. 

In line with these theoretical expectations, significant 
increases of (fc) with N and decreases of p with are 
apparent from Table I. As we will be comparing net- 
works generated from clusters of different size, it should 
be remembered that these basic features of the networks 
are changing. 



A. Small world properties 

One of the key properties when deciding whether a 
network behaves as a 'small world' is the average separa- 
tion between nodes, l^ve- More precisely, for each pair of 
nodes in the network the shortest pathway between them 
is found and then the average of the number of steps in 
these pathways is taken. Values of lave are given in Table 
I and its dependence on network size is depicted in Fig. 
2(a). Zavc seems to be growing sub- logarithmically with 
^minj although the growth is non-monotonic — such size 
effects are typical of small clusters. 

Two other quantities that reflect the properties of the 
shortest pathways between nodes are given in Table I. 
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FIG. 2: The dependence of (a) the average separation between 
nodes (in steps) and (b) the clustering coefficient, C, on the 
size of the network, nmin. The data points for LJ clusters with 
from 7 to 14 atoms are compared to the values expected for a 
lattice (Eq. (7)) and an Erdos-Renyi (ER) random graph (Eq. 
(8)), and to those for randomized (rand) networks with the 
same degree distribution as the clusters, as labelled. Values 
of both Ci and C2 are plotted for the LJ clusters. 



The diameter, d, is the number of steps in the longest 
such pathway, and the radius, r, is the number of steps 
in the longest pathway for the node for which this is a 
minimum, i.e. all the other nodes are within r steps of 
that node. The small values of Zave, r and d would suggest 
that the network is showing the 'small world' effect, but 
a more sophisticated analysis is also needed. 

The usual approach to decide whether the average sep- 
aration between nodes shows behaviour typical of a ran- 
dom graph or of a lattice is to test whether l^ve is a 
logarithmic (log n„) or power-law {nl^'^ where d is the di- 
mension of the lattice) function of the network size, «„, 
respectively. However, the situation is not that simple 
for our networks, because both the dimension of config- 
uration space and the average degree increases with the 
network size. Therefore, to test for lattice-like behaviour 
we need to compare the network, not to a lattice of fixed 
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dimension and increasing extent, but of fixed extent and 
increasing dimension. It is of fixed extent because tlie 
displacement of an atom that is required to reach tlie 
nearest permutational isomer of any minimum should al- 
ways be of the order of the atomic size. 

For a 3-/V-dimensional hypercubic lattice with L lattice 
points along each coordinate direction, the average sep- 
aration between lattice points is approximately the sum 
of the average difference for each of the 2>N coordinates. 
More specifically, 
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where the initial factor is to exclude the zero length paths 
between a node and itself from the average and 
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Hence l]^^^ = N{L~1){L+1)/L^'' Therefore, the average 
separation between lattice points scales linearly with the 
dimension of the system. The dependence of Z^v" '^^ the 
number of lattice points can be obtained by substituting 
for N in the above equation using niatt = , leading 
to 



,latt + 
'ave = ori^„ r logniatt- 



SLlogL 



(7) 



Thus, even in the case of a lattice, the average separation 
scales logarithmically, because of the increasingly high 
dimensionality of configuration space. 
For a random graph 
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log Uv 
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so for our networks, if this equation were obeyed we 
would expect a sublogarithmic scaling with size. In fact, 
if we substitute the result from Eq. (3) we get 



^ave — 



log n„ 



log (2a/ a) + log (log nmin) 



(9) 



The apparent sublogarithmic behaviour of Fig. 2(a) 
would suggest that as with the Watts-Stogatz small- 
world networks, the scaling of Zavc is similar to a ran- 
dom graph, but more conclusive evidence in favour of 
this conclusion comes if we compare Equations (7) and 
(8) to our data. To apply Eq. (7) a value for L was first 
obtained using the number of minima for LJ14, which 
gives L = 1.220. A much better fit to our data is obtained 
from the random-graph expression, confirming our initial 
assessment. 

This scaling behaviour may seem somewhat surprising, 
since the connections between minima are based on the 
adjacency of the associated basins in configuration space, 
which would perhaps initially suggest a more lattice-like 
picture of configuration space. Furthermore, there are 
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FIG. 3: Tlie cumulative distribution for tlie number of nodes 
tliat liave more than k connections. Tlic curves correspond to 
clusters of different sizes, as labelled. An additional straight 
line with slope —(7 — 1), where 7=2.78, has been plotted to 
emphasise the power law tail. In the inset the cumulative 
probability distribution for the 12-, 13- and 14-atom clusters 
is plotted against k/{k), to bring out the universal form of 
the distribution. 



no obvious equivalents of the random linkages between 
distant nodes that cause the small-world behaviour in 
the Watts-Strogatz networks. As we will see later, the 
origin lies elsewhere. 

Another feature of many networks is a strong degree 
of local structure, as measured by the clustering coeffi- 
cient. The clustering coefficient of node z, c^, is defined 
as the probability that a pair of neighbours of i are them- 
selves connected. However, when extending this concept 
to the whole graph there are two definitions in common 
usage. The first, C\, is simply the probability that any 
pair of nodes that have a common neighbour are them- 
selves connected, and the second, C2, is the average of 
the local clustering coefficient: C2 = X^iCi/A^min- The 
difference between these two definitions is the relative 
weight given to nodes with different degree. High degree 
nodes make a larger contribution to Ci because there will 
be more pairs of nodes that have a high-degree node as a 
common neighbour, whereas all nodes contribute equally 
to C2. Typically, C\ < C2 because, as is the case here, 
higher degree nodes tend to have lower values of Cj. 

Values of C'l and C2 are given in Table I and are com- 
pared to that for an Erdos-Renyi random graph (C^^ = 
p) in Figure 2(b). It is apparent that the networks for 
the larger clusters have significantly more local structure 
than the random graphs. For LJ14 Ci and C2 are 7.53 
and 35.90 times larger than C^^, respectively. At small 
cluster sizes C ~ Cer, mainly due to the large values of 
p — there is going to be little difference from a random 
graph when the network is almost fully connected — but 
as the networks become more sparse, they show increas- 
ing local correlations. This feature of the networks is 
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potential energy of minimum 

FIG. 4: (Colour online). The dependence of the degree of a 
node on the potential energy of the corresponding minimum 
for LJi4. The data points are for each individual minimum 
and the solid line is a binned average. 

unsurprising given that the connections are based on ad- 
jacency in configuration space. 



B. The degree distribution 

Another very important property of a network is its 
degree distribution. These distributions are illustrated 
for our networks in Figure 3. Remarkably, as the size 
of the clusters increase a clear power-law tail develops. 
Such networks are typically said to be scale-free, even 
though the power-law behaviour does not extend back to 
small values, and so the distribution is not totally invari- 
ant to rescaling.'*^ In fact, the initial flat portion of the 
cumulative degree distributions in Fig. 3 implies that the 
degree distribution itself has a maximum that lies near 
to, but below, (fc). Interestingly, the inset to Fig. 3 sug- 
gests that the distribution is tending to a universal form 
independent of cluster size. The value of the exponent 
in the power law is 2.78, which is quite similar to other 
scale- free networks (most are between 2 and 3). 

This result implies that the networks are very hetero- 
geneous, with a few hubs having a very large number 
of connections, but with most having something near to 
the average. It also explains the origin of the small-world 
behaviour, as scale-free networks are known to have very 
small values of Zavc- This is because the hubs provide 
natural conduits for short paths between nodes. In fact, 
the scaling of ^avc with network size can be even slower 
than for Erdos-Renyi random graphs. ^^'^^ 

The obvious next question to ask is what features dif- 
ferentiate those minima with large numbers of connec- 
tions from those with only a few. What is the source 
of the heterogeneity? Most of the models devised to ex- 
plain the scale-free properties of networks are dynamic 
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FIG. 5: (Colour online). (Cumulative) probability distribu- 
tions for the energy of the minima for LJ14. The distributions 
are compared to Gaussian forms (dotted lines). 



in nature. For example, in Barabasi and Albert's model 
it arises from the preferential attachment of new nodes 
to those with higher degree at each step in the growth 
process. However, our inherent structure networks are 
static in nature, and are just determined by the form of 
the interatomic potential and the number of atoms. 

One of the few scale-free models not involving network 
growth is the fitness model of Calderelli et ai, in which 
each node is assigned a fitness parameter that controls 
its properties. For a suitable combination of the distri- 
bution of the fitness parameter and the fitness-dependent 
linking probability, scale- free networks can result. 

The most obvious parameter that might control the 
properties of our nodes is the potential energies of the 
associated minima. In Figure 4 we show the dependence 
of the degree of a node on its potential energy for our 
largest network. There is a strong correlation; the global 
minimum is the most connected and the degree decays 
to lower values as the energy increases. This picture is 
typical of the data for other clusters, and only for LJg 
and LJio is the global minimum not the most highly- 
connected hub; in fact it is the fifth lowest-energy mini- 
mum in both these exceptions. The fcmax values in Table 
I indicate the highly-connected nature of the hubs. The 
biggest hub is always connected to more than 50% of the 
network, and the particularly high value for LJ14 (76%) 
explains the non-monotonic behaviour of /avo at = 14. 

The combination of k{E) illustrated in Fig. 4 with the 
probability distribution for the potential energy of the 
minima, p{E), must lead to a power-law tail for p(fc), 
and so this might be one avenue to explain the scale- 
free behaviour of our networks. Indeed, the form of 
p{E) has been extensively investigated for a variety of 
systems, and a physical motivation given for its form. 
Based on the minima sampled in the liquid state of bulk 
matter, the distribution of minima has been found to 
be a Gaussian function of the potential energy.^'^'-'^'* This 
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form can be justified from the central limit theorem if it 
is assumed the system can be divided into independent 
subsystems. 

Our data can also be well-approximated by a Gaus- 
sian (Figure 5), however there are deviations in the tail 
regions. In particular, p{E) decays more slowly than a 
Gaussian at low energy. This deviation is not so surpris- 
ing, since the analysis is expected to be most applicable 
where there is a quasicontinuous distribution of minima, 
rather than for the more discrete spectrum of minima 
in the low-energy tail that is associated with the solid 
form of the clusters. However, this part of the distribu- 
tion is particularly significant as it determines the prop- 
erties of the highly-connected minima in the scale-free 
tail. Indeed, a more detailed analysis shows that the 
non-Gaussian nature of the low-energy tail of p{E) com- 
bines with the faster than exponential increase of k{E) 
to produce the power-law tail for p{k). However, it is no 
more clear why these functions should have these precise 
forms than why p{k) has a power-law tail, so relatively 
little new physical insight has been obtained. 

Probably, a more fruitful way of understanding the 
heterogeneous nature of the networks is by thinking in 
terms of the hyperareas of the basins of attraction sur- 
rounding a minimum. The high connectivity of the low- 
energy minima is likely to stem from their large basins 
of attraction, because this gives them long basin bound- 
aries, and hence the likelihood of a greater number of 
transition states on these boundaries. Indeed, there is 
evidence that for LJ38 the basins of attraction rapidly 
decrease in area with increasing potential energy.^'^ Some 
statistics are also available for the probability of reaching 
the global minima of L J clusters from a random starting 
point, which provide a direct measure of the relative area 
of the global minimum. ^^-^^ These studies indicate that 
this probability is significantly larger than l/nmin, and 
thus that they have large basins of attraction. We intend 
to pursue this idea further in future work, in particular to 
examine in detail whether the relationship between basin 
area and connectivity is as we have suggested. 

These ideas lead to a hierarchical picture of the en- 
ergy landscape, where large basins of attraction are sur- 
rounded by smaller basins, which in turn are surrounded 
by smaller basins, .... A useful way to envisage this 
heterogeneity of basin areas is by analogy to an Apol- 
lonian packing, a two-dimensional example of which is 
shown in Fig. 6. The initial configuration for this pack- 
ing is three mutually touching discs. In the interstice 
between these three discs a new disc is added that itself 
touches the three discs (the so-called inner Soddy circle). 
In the three new interstices created, three more discs are 
added. This procedure is continued iteratively to gener- 
ate a space-filling packing of discs that is a fractal^^ of 
dimension 1.3057.'"'^ 

As with the basins of attraction that divide configu- 
ration space, the discs in the Apollonian packing have a 
very heterogeneous distribution of areas that is organized 
hierarchically. In fact, the probability distribution of the 




FIG. 6: An Apollonian packing of disks. 



areas of the disks is a power-law. Furthermore, if we 
consider the Apollonian packing as a network where each 
disc corresponds to a node and nodes are connected by an 
edge if the corresponding discs are in contact, then this 
network has a scale-free degree distribution.^*'^^ More- 
over, many of the topological properties of this "Apol- 
lonian network" are similar to the inherent structure 
networks.^* 

This analogy helps us to understand why our initial 

intuition namely, that the inherent structure networks 
would be lattice-like, because they are based on the adja- 
cency of basins — was wrong. However, there are of course 
limitations to the analogy. In the Apollonian packing 
space is filled by an infinite number of discs, whereas 
configuration space is divided up into a finite number of 
basins surrounding the minima. Therefore, it's proba- 
bly better to compare the ISN to an Apollonian packing 
where the iterative addition of discs has only been ap- 
plied a finite number of times. ^® 

Having outlined the basic properties of the inherent 
structure networks, it is useful to compare our results to 
the other studies looking at the connectivity of configura- 
tion space. The scaling behaviour of Zave has been char- 
acterized in both Scala et al.'s study of lattice polymers'*^ 
and Wuchty's analysis of RNA.*"^ In Ref. 42 networks of 
different size were generated by considering both poly- 
mers of different length and with different end-to-end 
distances. Based on a comparison of Zave to that for a 
two-dimensional lattice and a random graph, they con- 
cluded that there is a crossover to logarithmic scaling for 
networks of size larger than 100. However, this is the 
wrong comparison, because the dimension of configura- 
tion space is of course, 2iV, where N is the number of 
monomers in the chain. Instead, the deviation from a 
power-law scaling is probably due to the effective higher 
dimensionality of conformation space for the longer poly- 



8 



mors. Consistent with this interpretation the degree dis- 
tribution shows no sign of heterogeneities, but is instead 
a Gaussian, because of the local nature of the move set. 
Each conformation occupies a similar volume in confor- 
mation space. 

Wuchty considers only one RNA sequence of a fixed 
length, but generates networks of different sizes through 
only including states that are within a certain variable 
energy of the groimd-state. In contrast to our graphs, 
the networks are therefore incomplete, /ave 

deviates sig- 
nificantly from that for a random graph, and the degree 
distribution has a exponential tail. 

We would suggest that both these networks are essen- 
tially behaving as lattices of high dimensionality, because 
the local move-sets only allow connections to conforma- 
tions that are nearby in configuration space, and there is 
no equivalent to the inherent structure mapping to group 
conformations into states with widely-differing basin ar- 
eas. 

By contrast, in their study of a protein network Rao 
and Caflisch found behaviour much more similar to our 
inherent structure networks with a degree distribution 
that had a power-law tail and with the low-lying free- 
energy minima having greater connectivity.^* This re- 
flects the conceptual similarities in the network repre- 
sentations noted in Section II and suggests the potential 
generality of the results we have found here. 



C. Randomization 

Before wc look at more detailed properties of the net- 
works, it is useful to have an appropriate random model 
to which to compare. For example, Dorogovtsev and 
Mendes have suggested that the values of C /C™- for our 
inherent structure networks simply reflect the degree dis- 
tribution and not any additional local structuring.^*^ 

Analytic expressions are available for some network 
properties of a random ensemble of networks with a given 
degree distribution. However, these do not apply when 
multiple edges and self-connections are excluded. In- 
stead, we used the switching algorithm to generate the 
ensembles of random networks. At each step in this 
procedure two edges are picked at random, A — B and 
C — D say, and then rewired to form two new edges, ei- 
ther A — C and B — D, or A — D and B — C. It is evident 
that this rewiring conserves the degree distribution. It is 
also straightforward to keep the constraint that there are 
no multiple edges or self-connections; rewiring steps that 
would break this constraint are simply not accepted. The 
procedure is started from the original network, and, after 
a suitable equilibration period, statistics were recorded. 
At each step the probability of each edge occurring is 
followed. From this probability distribution a variety of 
network properties can be calculated. However, for other 
properties, e.g. Zavei that have to be explicitly calculated 
for a specific network, averages over 25 independent (as 
judged by the autocorrelation function of the clustering 



TABLE II: Some properties of the randomized networks pro- 
duced from the inherent structure networks by the switching 
algorithm.^°'^i 



N 


^ave 


Ci 
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1.479 


0.625 


-0.408 


-0.299 


9 


1.706 


0.492 


-0.304 


-0.125 


10 


2.043 


0.332 


-0.158 


-0.088 


11 


2.071 


0.219 


-0.108 


-0.063 


12 


2.247 


0.124 


-0.115 


-0.085 


13 


2.338 


0.071 


-0.106 


-0.109 


14 


2.302 


0.032 


-0.068 


-0.119 



coefficient) random networks were taken. 

Some properties of these randomized networks are 
given in Table II. For example. Figure 2(b) shows that 
although the randomized networks have a significantly 
higher clustering coefficient than the Erdos-Renyi ran- 
dom graphs, the values for the inherent structure net- 
works are yet still higher. This result disproves Doro- 
govtsev and Mendes' assertion^" and shows that there 
are additional sources of local structuring, over and above 
that due to the degree distribution. 

By contrast, in Figure 2(a), which compares /avc for the 
cluster networks and their randomized versions, only very 
small differences are apparent, showing that this property 
primarily reflects the degree distribution of the networks. 
However, it is noticeable that ^av"^ is always systemati- 
cally smaller, albeit by a small amount, except for LJg. 
This reflects the general correlation between higher clus- 
tering and longer path lengths, as seen, for example, in 
the extremes of the Erdos-Renyi random graphs (low C, 
small Zave) and lattices (high C, large Zave)- In a net- 
work with higher clustering, edges are 'wasted' connect- 
ing nodes that arc already close, which might have been 
used to connect distant nodes and hence lower ?ave more 
significantly. 

In the following we will only show data for the ran- 
domized networks when they show significant differences 
from the inherent structure networks. Typically, there 
will only be very small differences for properties related 
to the shortest paths on the networks, as with ^avei but 
more significant differences appear for properties associ- 
ated with clustering and network correlations. 

D. Degree dependent properties for LJ14 

The extremely heterogeneous nature of our networks is 
reflected in their scale-free degree distribution. To probe 
this heterogeneity further, it is also helpful to look at 
how properties vary from node to node, and particularly 
how they depend on the degree of a node. We will do 
this for our largest network, LJ14, but results for smaller 
networks are similar. Because of the large size of this 
network, there is often considerable scatter in the data. 
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FIG. 7: (Colour online). The dependence of the local cluster- 
ing coefficient of a node on the degree. The data points are for 
each LJi4 minimum, and the the sohd hues are the average 
values for a given k. For comparison, the function and 
c(fc) for the randomized networks and an Erdos-Renyi (ER) 
random graph have also been plotted. 



FIG. 8: (Colour online). The dependence of the average sep- 
aration of a node in the network from all the others on the 
degree. Each LJ14 minimum is represented by a data point 
and the solid line is the average value for a given k. I gmin is 
the average separation between the global minimum and the 
rest of the nodes. 



To make the trends more clear, averages of the properties, 
e.g. over nodes with the same degree, are also usually 
presented. 

The dependence of the local clustering coefficient on 
degree is shown in Figure 7. Like many other net- 
works their is a clear association between higher de- 
gree and lower clustering. ^''"^^"^^ This property has of- 
ten been attributed to a hierarchical structuring of 
the network. ^^'^"^ Notably many of the deterministic 
scale-free networks, ^'^■^^'^^ including the Apollonian net- 
work mentioned in Section IIIB,^^ follow the power law 
Ci{ki) ~ . At large k the LJ14 network roughly fol- 
lows this law, but at low k the Ci values do not increase 
as fast as expected by this law. 

It is noteworthy that Ci{k) for the randomized net- 
works has a similar type of the dependence on fc, as for 
the LJ clusters. This is somewhat surprising since for a 
random uncorrelated graph in which multiple edges and 
self-connections are allowed, the local clustering coeffi- 
cient is independent of degree, no matter what the de- 
gree distribution.^^ The behaviour of c-™^(fc) is due to 
the correlations induced by the absence of multiple edges 
and self-connections^"'^^ and is something we will return 
to in the next section. 

At small k the local clustering coefficient for the ISN 
is roughly 50-100% larger than for the randomized net- 
works. This extra local structuring most likely reflects 
the spatial embedding of our networks. Low degree nodes 
have small basins of attraction, and so are only connected 
to basins that are close in configuration space. This spa- 
tial localization of the connections leads to high cluster- 
ing. By contrast, the large degree nodes have large basins 
of attraction, and so the minima surrounding them can 
be distant in configuration space, and are not especially 
likely to also be connected to each other. Indeed, the 



clustering coefficient for the global minimum is similar 
to the value expected for an Erdos-Renyi random graph. 

In Figure 8 we show U, the average separation of a 
node i from the rest of the network. It is clear that the 
global minimum plays a central role in the network. It is 
the node that is on average closest to all the other nodes 
in the network. Similarly, for all our networks the global 
minimum is one of the nodes for which all other nodes 
are within the radius of the network (Table I). By con- 
trast, low degree nodes are more on the periphery of the 
network. It is also apparent from Fig. 8 that the minima 
separate themselves into sets depending on how far they 
are away from the global minimum. For those that are 
directly connected to the global minimum, ^gmin + 1 is an 
upper bound to U and U would take that value if all the 
shortest paths involving that node passed through the 
global minimum. Similarly, the upper bound for those 
two steps away from the global minimum is ^gmin + 2, 
and so on. In fact there are only a few minima that are 
three steps away from the LJ14 global minimum. 

The vertex betweenness measures the fraction of the 
shortest paths between all nodes in the network that pass 
through a particular vertex. Similarly the edge between- 
ness measures the fraction of these paths that pass along 
a particular edge in the network. These quantities pro- 
vide a measure of the importance of a node or edge to a 
network, particularly for dynamical processes that occur 
on the network, such as the passage of information. For 
scale-free networks, there is typically a strong correlation 
between the vertex betweenness and the degree, and 
Figure 9(a) shows that the inherent structure networks 
show a similar behaviour. The vertex betweenness has 
roughly a power-law dependence on k, and so given the 
scale-free degree distribution, this also implies that the 
probability distribution for the betweenness also has a 
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FIG. 9: (Colour online). The dependence of (a) the vertex 
betweenness on the degree and (b) the edge betweenness on 
the product of degrees for LJ14. Each data point represents 
(a) a minimum or (b) an edge, and the solid line the average 
value for a given (a) k or (b) range of kikj. 



power-law tail.'''^ Again the global minimum plays a cen- 
tral role in the network. For LJ14 59.8% of the shortest 
paths pass through this node. 

For scale-free networks, there is also usually a corre- 
lation between the edge-betweenness and some measure 
of the degrees at either end of the edge, but it is much 
weaker than for the vertex betweenness.''^ In Figure 9(b) 
we see such a correlation with the product of the degrees 
at either end of the edges for the LJ14 network. However, 
the correlation is fairly weak and there is much more scat- 
ter in the data. The edge with maximum betweenness 
connects the two highest-degree nodes and 0.22% of the 
shortest paths pass along it. 



E. Network Correlations 

To further understand the structure of our networks, 
it is important to go beyond just local properties of a 
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FIG. 10: (Colour online). The dependence of the average 
degree of the neighbours of a node on the degree of that node. 
Each LJ14 minimum is represented by a data point and the 
solid line labelled 'LJ' is the average value for a given k. Also 
included is knn{k) for the randomized networks, a line with 
the exponent expected for an Apollonian network,^* and a line 
at k^^°'''' , the value expected for an uncorrelated network. 



node to study the correlations between the properties 
of the nodes. Indeed correlations can have an impor- 
tant influence on other network properties, and have 
a strong effect on dynamical processes that occur on a 
network. ^''"^^ 

Most commonly, correlations with respect to degree are 
investigated. For example, in Figure 10 we show how fc„„, 
the average degree of the neighbours a node, depends on 
the degree of the node. If the network were uncorre- 
lated, one would expect fc„„ to be independent of degree 
and to take a value k"^"" = However, knn{k) 

has a negative slope showing that our networks are dis- 
assortative by degree, i.e. high degrees are more likely 
to be connected to low degree nodes than expected, and 
vice versa. Conversely, in assortative networks nodes are 
more likely to be connected to those with similar degrees. 
Interestingly, fc„„(fc) has an approximate power-law de- 
pendence on fc, as seen for other scale-free networks, "^^ 
with an exponent that is very similar to that of an Apol- 
lonian network.^* 

Another way to measure such correlations is through 
the assortativity coefficient introduced by Newman. ^^■'^^ 
It is a two-point correlation function of the properties at 
either end of an edge and is usually normalized by the 
value expected for a perfectly assortative network. It is 
defined as 



St), - {s),{t). 



i^i) e,as 



(^)e(i). 



(10) 



where s and t correspond to the property of interest at ei- 
ther end of an edge, e denotes that the averages are over 
all edges and assort that the average is for a perfectly 
assortative network. A positive value is expected for an 
assortative network (with an upper bound of 1), zero for 
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FIG. 11: (Colour online). The dependence of the potential 
energy of the neighbours of a minima on its energy. The 
data points are for every LJ14 minimum, and the lines binned 
averages for both the LJ14 network, and randomized versions, 
as labelled. 

an uncorrclatcd network and a negative value for a disas- 
sortative network. Negative values have been found for 
technological networks, such as the internet and world- 
wide web, and biochemical networks, such as the network 
of protein-protein interactions and metabolic networks, 
whereas strongly positive values have been found for so- 
cial networks. ^^''^^ Values of are given in Table I and, 
as with knn {k) indicate the disassortative nature of our 
networks. The values found for the larger clusters are 
quite similar to those for technological networks. '^^''^^ 

The origin of this disassortativity has been explored 
in most depth for the internet. In particular, it 
was found that the exclusion of multiple edges and self- 
connections can be a significant source of disassortativity. 
Indeed, when we compare fc„„ and to the values for 
the randomized networks (Fig. 10 and Table II), we find 
very similar levels of disassortativity. 

In an uncorrelated network, the probability that a ran- 
domly chosen edge will connect two particular nodes is 
/ci/c2/2?ig(jggs and, hence, the total number of expected 
connections between the two nodes is fcifc2/2nedgos- For 
the LJ14 network this gives an expected value of 33.5 
edges between the two highest degree nodes. As only 
one is possible, the network appears disassortative. 

Correlations also have a significant effect on the 
clustering. As indicated by the expressions above, 
the probability that the neighbours of a node are con- 
nected depends sensitively on the degrees of the neigh- 
bours, with high degree nodes being much more likely 
to be connected. Hence the local clustering coefficient is 
higher for nodes with larger fc„„ , and so even for the ran- 
domized graphs the local clustering coefficient decreases 
as k increases (Fig. 7). 

We have also explored correlations between the poten- 
tial energies of connected minima. Values of r£; are given 
for both the inherent structure networks and the random- 



ized versions in Tables I and II, respectively, and the be- 
haviour of Enn , the average energy of the the minima to 
which a minimum is directly connected, is shown in Fig- 
ure 11. Because of the correlations between k and E (Fig. 
4), one might also expect the networks to be disassorta- 
tive with respect to the energy of the minima. However, 
for N > 9 positive values of are found, whereas the 
values for the randomized networks are always negative, 
and of roughly similar mag nitude to rf""^. This differ- 
ence shows that the potential energy of the minima plays 
an important role in the structure of the network. 

At low energies, Enn{E) has a gentle negative slope 
similar to that for the randomized graphs, and probably 
mainly reflects the exclusion of self-connections and mul- 
tiple edges. The low-enegy hubs have to be connected 
to higher-energy minima because of the constraint that 
all the neighbours must be different. For example, the 
LJ14 global minimum is connected to 3201 of the 4196 
nodes. Even if these were the 3201 lowest-energy nodes 
(i.e. maximally assortative given the exclusion of multi- 
ple edges), the average energy of its neighbours would 
only be a little lower (0.28e) than the actual value. How- 
ever, for the higher-energy nodes En„{E) has a strongly 
positive slope, because of the preference for linking to 
minima with similar energy. 

Although uncorrelated energy landscapes have been 
studied thcoretically,^^'^^ the assumption that the energy 
of connected states is independent is usually unrealistic 
and gives rise to landscapes that are extremely rough. 
(Indeed the first step in incorporating greater realism in 
these models is often to include correlations. ^'^^^^) For 
landscapes associated with the configurations of atomic 
systems, the landscapes are much smoother, simply be- 
cause connected states are likely to be structurally similar 
and hence have similar energies. However, it should be 
remembered that the measures of correlations we have 
used here are quite local, and apply to just the immedi- 
ate neighbourhood of a basin. When energy landscapes 
are inferred to be rough, e.g. because the system is a 
good glass-former and easily gets stuck in traps on the 
landscape,'* this is usually referring to a lack of correla- 
tions at larger length scales. 



IV. DYNAMICAL IMPLICATIONS 

So far, we have made a detailed characterization of 
the network of connections between the minima on a po- 
tential energy surface, and these results provide impor- 
tant insights into the fundamental organization of energy 
landscapes. But it is also important to understand the 
implications of this topology for the dynamics of a sys- 
tem. 

Clearly, the small-world nature of the inherent struc- 
ture networks has important consequences for the the- 
oretical understanding of searching configuration space. 
For example, the Levinthal paradox in protein folding 
highlights the huge number of possible configurations of 
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a protein and the seeming impossibility of reaching the 
native state in a random search through these configura- 
tions, in contrast to the protein's actual ability to fold.® 
Similar paradoxes can of course be formulated for any 
system with a large configuration space, be it a macro- 
molecule, cluster or bulk liquid, that is able to routinely 
find a preferred structure. More formally, the exponen- 
tial increase in the size of the search space underlies the 
classification of most problems involving the global op- 
timization of a configuration as "NP-hard";^^ i.e. there 
is no known algorithm that is guaranteed to solve the 
problem in polynomial time. 

Our results offer part of the answer to these paradoxes. 
Even for high-dimensional configuration spaces with huge 
numbers of accessible states, the number of steps in a 
pathway from a random configuration to the target state 
remains modest due to the favourable scaling with size 
(Eq. (9)). Of course, finding that particular pathway 
in the absence of additional guiding information can be 
extremely hard, although there is some evidence that a 
local knowledge of the topology can be of some help.^®'^^ 
Indeed, the emphasis of most proposed solutions to the 
Levinthal paradox is on those features of the energy land- 
scape, such as funnels, that guide the system in the 
right direction in its descent down the energy landscape, 
and whose presence differentiates those systems that are 
able to find the target state, from those that get stuck in 
the morass of possible configurations. 

The dramatic consequences of this small-world be- 
haviour are evidenced in examples where the energy land- 
scape has a favourable topography. For example, for LJ55 
it has been estimated that there are 10^^ minima. How- 
ever, the basin-hopping global optimization algorithm is 
able to find the global minimum from a random starting 
point after sampling on average only approximately 150 
minima. 

The scale-free nature of the inherent structure network 
will also have important implications for a system's dy- 
namics. Interestingly, at the centre of Leopold et aPs 
original definition of a folding funnel was the concept of 
a "convergence of multiple pathways" at the native state 
of a protein.^^ Interpreting this idea in network terms, 
it seems to be suggesting that the native state would be 
a hub in the network. For our inherent structure net- 
works we see such convergence with the global minimum 
connected to a significant fraction of the other minima, 
and similar results have been obtained in Rao et aZ.'s in- 
vestigation of the connectivity of configuration space for 
polypeptide chains.'^'' Of course, there is also a strong to- 
pographical component inherent in the idea of a funnel, 
namely that as the system goes downhill it is becoming 
closer to the native state, and it is this feature that has 
received the most emphasis more recently, perhaps to the 
neglect of the potential focussing aspects of the topology. 

A wide variety of dynamical processes have been stud- 
ied on complex networks. The dynamics of interest here 
is the molecular dynamics on the potential energy sur- 
face, as dictated by the forces on the system. The most 



relevant work to this from the networks literature is that 
examining diffusion on scale- free networks,*^ where it was 
found that this topology led to more favourable scaling of 
the transit time between nodes, particularly when these 
nodes had high-degree. 

The effect of the topology on the dynamics on an en- 
ergy landscape can be illustrated for a simple model land- 
scape where topographical features are eliminated. The 
model consists of a a flat landscape where all the minima 
are identical, except for their connectivity, i.e. they have 
the same energies and vibrational (and hence thermo- 
dynamical) properties, and furthermore, all the transi- 
tion rates between connected minima take the same value 
The average residence time in a minimum i is then 
1/kikK As the equilibrium probability of being in a min- 
imum i is simply l/rimin, the frequency with which this 
minimum is visited is A;,fcV'^inin- The first passage time 
for encountering a minimum therefore decreases with de- 
gree. The topology of the inherent structure network has 
a focussing effect directing the system more rapidly to 
the highly-connected hubs. The effects of the topology 
on the dynamics will be explored in more detail in future 
work. 



V. CONCLUSION 

In this paper we have provided a detailed characteri- 
zation of the inherent structure networks for a series of 
small Lennard-Jones clusters. These networks show the 
mixture of local order (as measured by clustering coeffi- 
cients) and small average separations between nodes that 
is characteristic of small-world networks. Furthermore 
as the size of the clusters increase, these inherent struc- 
tTirc networks develop a clear power- law tail to the degree 
distribution and have a whole variety of properties that 
are typical of scale-free networks. 

However, in contrast to most scale-free networks the 
origin of this scale-free topology is not due to some form 
of preferential attachment during network growth.^''' In- 
stead, the network heterogeneity reflects the topography 
of the energy landscape with the degree of a node strongly 
correlated to the potential energy of the associated min- 
imum. This correlation most likely arises because the 
low-energy minima have larger basins of attraction, and 
so are likely to have more transition states located on 
their basin boundaries. Apollonian networks provide a 
model of how such basins could be organized to give rise 
to a scale-free network. Hence, the discovery of the scale- 
free character of inherent structure networks may have 
profound implications for our understanding of the way 
the inherent structure mapping divides up configuration 
space, and raises the possibility that configuration space 
is tiled by a fractal, Apollonian-like packing of basins. 

The results presented here for the inherent structure 
networks are for systems interacting with a particular 
potential and of small size, so it is natural to ask how 
general are the results. Firstly, we see no obvious rea- 
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son why the Lcnnard- Jones potential should be "special" , 
producing scale- free networks, while the landscapes as- 
sociated with other potentials have different topologies. 
However, it is conceivable that differences might arise as 
a result of the orientational degrees of freedom associated 
with molecular (not atomic) systems or the constraints 
of chain connectivity associated with polymers, and we 
intend to explore this further in future work. The work 
of Rao and Caflisch is also of interest in this regard, be- 
cause the networks they constructed to represent the con- 
figuration space connectivity of polypeptides also showed 
scale-free behaviour. 

Secondly, because of the exponential increase in the 
number of minima and transition states with system size, 



the analysis we presented here is inevitably limited to 
relatively small sizes. Furthermore, there is currently no 
known method to construct a statistical representation 
of the inherent structure network from an incomplete 
sampling of the network. One potential way to anal- 
yse systems of larger size is to use a more coarse-grained 
division of configuration space than provided by the in- 
herent structure mapping, such as that used by Rao and 
Caflisch.'*^ Another more indirect approach would be to 
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